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Abstract 

We propose a novel efficient algoritlim to solve Poisson equation in irregular two dimensional 
domains for electrostatics. It can handle Dirichlet, Neumann or mixed boundary problems in 
which the filling media can be homogeneous or inhomogeneous. 

The basic idea of the new method is solve the problem in three steps: (i) First solve the 
equation V • D = p. The inverse of the divergence operator in a restricted subspace is found to 
yield the electric flux density D by a fast direct solver in 0{N) operations. The D so obtained 
is nonunique with indeterminate divergence- free component. Then the electric field is found by 
E = D/e. But V X E = for electrostatic field; hence, E is curl free and orthogonal to the 
divergence free space, (ii) An orthogonalization process is used to purify the electric field making 
it curl free and unique, (iii) Then the potential (j) is obtained by solving Vcj) = — E or finding the 
inverse of the gradient operator in a restricted subspace by a similar fast direct solver in 0{N) 
operations. 

Treatments for both Dirichlet and Neumann boundary conditions are addressed. Finally, the 
validation and efficiency are illustrated by several numerical examples. Through these simulations, 
it is observed that the computational complexity of our proposed method almost scales as 0(A), 
where A is the triangle patch number of meshes. Consequently, this new algorithm is a feasible 
fast Poisson solver. 
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1 Introduction 



The Poisson equation occurs in the analysis and modehng of many scientific and engineering problems. 
In electrostatics, Poisson equation arises when finding the electrostatic potential of an electric field 
in a region with continuously distributed charges [1 . It is often solved in micro- and nano-electronic 
device physics [2], as well as in electronic transport and electrochemistry problems in terms of the 
Poisson-Boltzmann equation [3 . In fiuid dynamics, Poisson equations are solved to find the velocity 
potential in a steady-state potential fiow of an incompressible fiuid with internal sources or sinks j4j. 
Moreover, Poisson equation is also encountered in finding the steady-state temperature in an isotropic 
body with internal sources [5]. 

An accurate and efficient solution of Poisson equation is critical in various areas. For example, 
in design optimization of nanodevices where quantum effects are significant, a widely used scheme is 
to solve the coupled Schrodinger-Poisson system self-consistently [H [7] , in which Poisson equation is 
solved repeatedly, and concomitantly with the Schrodinger's equation. As such, the computational 
load for solving Poisson equation is always of concern. 

There are two main classes of solvers for linear systems from Poisson equation: direct and iterative. 
One of the direct solvers is the fast Poisson solver based on fast Fourier transform (FFT) [8 . Indeed, 
this method is extremely efficient when the solution regions are simple and regular geometries with 
regular grids, such as rectangular regions, 2-D polar and spherical geometries [9], and spherical shells 
[To]. Since practical problems usually involve complex geometries, there have been many research 
works on seeking alternative methods. The multifrontal method with nested dissection ordering [11] 
is the most efficient direct method that can deal with complex geometries. Its key idea relies on 
partitioning the domain using a nested hierarchical structure and generating the LU decomposition 
from bottom up to minimize fill-ins. Typically, the computational complexity of the multifrontal 
method is of 0{N^'^) in two dimensions where N is the dimension of the matrix. 

The other class of solvers, the iterative ones, are more favorable when large systems are solved. A 
popular one is the approach based on integral equation techniques and accelerated by fast multipole 
method (FMM) |12l[T3l[14l [15] . It can achieve 0{N) complexity when the underlying Green's function 
is available and amenable to factorization. Nevertheless, the underlying Green's function is difficult 
to obtain unless involved problems have constant coefficient in free-space or separate geometries with 
simple boundary conditions. 

Multigrid method is one of most effective preconditioning strategies for iterative Poisson solvers 
[T6[ IT7[ lis] . Although this method takes advantage of fine meshes and coarse meshes and can achieve 
nearly optimal efficiency in theory, it is difficult to implement in a robust fashion because it demands 
a set of hierarchical grids of different density, which is not convenient in many real world problems. 

In this paper, we present a new efficient numerical solution of Poisson equation for arbitrary two- 
dimensional domain with homogeneous or inhomogeneous media. Unlike traditional Poisson solvers, 
where the potential is found directly, this new method solves for electric ffux density D first. The 
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electric flu4j is first approximated by the RWG basis. It uses the loop- tree decomposition of the RWG 
approximated vectorial electric flux which is a quasi-Helmholtz decomposition method developed in 
computational electromagnetics (GEM). With this technique, the electric flux is expressed as the 
combination of the loop space (subspace) (solenoidal or divergence free) part and the tree space 
(subspace) (quasi-irrotational) part. These two spaces, however, are non-orthogonal to each other, 
unlike the case of a pure Helmholtz decomposition. 

Expanding the electric flux with two sets of basis functions: loop and tree basis functions, we can 
solve for the electric flux by a two-stage process: First, to find the tree-space part, a matrix system 
is derived from the differential equation, V • D = p, and then solved by a fast tree direct solver with 
0{Nt) complexity, where Nt is the total number of tree basis functions. The obtained electric flux is 
nonunique as a divergence free component can always be added to it. 

Second, the electric field E = D/e is obtained but it is incorrect because E should be curl- free 
but the so-obtained electric field is not curl-free. It is polluted by the loop-space components. The 
loop-space part of E is purged by a projection procedure, which is iterative. Having acquired the 
curl-free electric field distribution, which is unique, we can readily get the potential distribution by 
solving E = — by the fast tree solver as well. 

In addition, we address the special treatments for Dirichlet and Neumann boundary conditions, 
respectively, which guarantee that the obtained solution is unique. Through numerical examples, we 
validate the feasibility and efficiency of the new method whose computational complexity almost scales 
as 0{N), where TV is the triangle patch number of meshes. 

To the best of our knowledge, this is the first time that quasi-Helmholtz decomposition techniques, 
e.g., loop-tree decomposition, are introduced to solve Poisson equation. This new method provides 
a feasible alternative for existing fast Poisson solvers. Moreover, the work about Poisson equation 
with Neumann boundary condition in homogeneous media has been reported in fi9\. Here, not only 
does this paper extend the method to inhomogeneous medium cases, but it also addresses the special 
treatment for the Dirichlet boundary condition. 

The organization of this paper is as follows. In Section [21 we introduce the basic problem definition 
and several relevant preliminaries. In Section [H we present the algorithm for solving Poisson problem 
with Neumann boundary condition. Next, we offer a special treatment for Dirichlet boundary condi- 
tion in Section HI Finally, in Section O we will verify the validation and illustrate the efficiency of the 
new method by several numerical methods. Gonclusions will be drawn in Section [6l 
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Figure 1: Schema for a typical Poisson problem. 

2 Preliminaries 

2.1 Two dimensional Poisson problems 

Assume inhomogeneous dielectric materials occupying a two dimensional bounded and simply con- 
nected region, with boundary T and normal n that points to the solution region as shown in Fig.[TJ 
The boundary T is composed of two parts: the first one, denoted by Tjj^ is imposed by Dirichlet 
boundary condition and the other part. Fat, is imposed by Neumann boundary condition. 

In this paper, we are interested in solving the 2-D Poisson equation arising from electrostatic 
problem: 

V • e^(r)V^(r) = -p(r)/eo for r G 1^ 

(/)(r) =(/)o(r) forrGFz, (1) 

|^(/>(r) = ^(r) for r G F^ 

where is the potential and p describes the charge density, and eo are the relative permittivity 
and free space permittivity, respectively. Suppose that the Dirichlet boundary consists of finite M 
distinct boundaries, Tjj = Ui^i then a fixed potential is prescribed on the boundary F^'^ 
for i = 1,2,...,M. Moreover, to complete the description of a well-posed problem, the Neumann 
boundary data ^(r) must be a square integrable function over the corresponding boundary [20j. 

Apparently, when Fat = 0, the above equation shrinks to a Dirichlet problem, which has unique 
solution. Similarly, it becomes a Neumann problem that is uniquely solvable (up to a constant) when 
Fi) =0. 

2.2 RWG basis 

The Rao- Wilton- Glisson (RWG) function [ST is the most popular choice of expansion function (also 
called basis function) in the computational electromagnetics (CEM) community. In this paper, we use 

-•^The "electric flux density" is also called the electric displacement fleld. We will call it "electric flux" for short. 
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Figure 2: The geometry description of a Rao- Wilton- Glisson (RWG) function. 



a normalized version by removing the edge length from the original definition. As shown in Fig. [21 
the normalized RWG function straddles two adjacent triangles, and hence, its description requires two 
contiguous triangles. The expression for the expansion function describing the field on the triangle is 

±^(r-r±) ifrGT± 
A„(r) = <; " (2) 

otherwise 

where ± denote the respective triangles, and are the vertex points and areas of the respective 
triangles, and are the supports of the respective triangles. 

Apparently, RWG function is the lowest order divergence conforming function which is defined on 
a pair of triangles. Such basis functions can be used to expand the electric fiux 

N 

D(r) = ^A„(r) (3) 

where N is the total number of normalized RWG basis functions. Here, D is defined by 

D(r) = -e(r)V^(r) (4) 

where 

e(r) = eoe^(r) (5) 

is the electric permittivity. 

The divergence of A^, which is proportional to the charge density associated the basis function, 
reads 



V-An(r) 



(6) 



otherwise. 

The charge density is thus constant in each triangle. The total charge associated the triangle pair, T+ 
and T~, is zero. This implies charge neutrality physically. Hence, the divergence conforming property 
of RWG function makes its representation of the field physical. 

2.3 Helmholtz theorem and loop-tree decomposition 

The well-known Helmholtz theorem [22] states that a vector field can be split into the form 

f (r) = V(/:) + V X V. (7) 
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The first term V^p is the irrotational (curl- free) part, and the second term V x v is the solenoidal 
(divergence- free) part. Low frequency problem has been an intense research area for the last decade in 
the CEM community. When the frequency is low, the current naturally decomposes into a solenoidal 
(divergence- free) part and an irrotational (curl- free) part. These two parts are imbalanced when the 
frequency becomes low. Hence, there is a severe low frequency breakdwown numerical problem when 
solving integral equations in which RWG function is normally used. One remarkable remedy is the 
well known loop-tree decomposition [23l [231 [25l [26l [271 [2H] , in which RWG basis is decomposed into the 
loop basis functions that have zero divergence, and the tree basis (or star basis) functions that have 
nonzero divergence. This is a quasi-Helmholtz decomposition because the tree expansion functions are 
not curl free. For short, we call the subspace spanned by RWG basis functions the RWG space, the 
subspace spanned by the tree basis functions the tree space, and that spanned by loop basis functions 
the loop space. 

Quite similar to the case of CEM, we can expand the electric flux 



where D/ and are the loop-space part and the tree space part respectively, (r) is a loop expansion 
function such that V • L^(r) = 0, and Ti(r) is a tree expansion function such that V • T^(r) 7^ 0, the 
numbers of the loop basis functions and the tree basis functions are Ni and Nt^ respectively. 

2.4 Loop basis and tree basis 



Figure 3: Linear pyramid basis function and loop basis function. Left: pyramid function. Right: loop 
basis function. 

A loop basis function is described by the surface curl of a vector function, namely. 



where the scalar function A(r), also referred as "solenoidal potential" [29 , is the linear Lagrange or 
nodal interpolating basis. As shown in the left part of Fig.[3l Ai(r) is a piecewise linear function with 
support on the triangles that have a vertex at the ith node of the mesh, attaining a unit value at 
node i, and linearly approaching zero on all neighboring nodes. This interpolating scalar functions is 
also intuitively called pyramid basis functions. Moreover, the right one of Fig. [3] illustrates the loop 




(8) 




L(r) = Vs X z A(r) 



(9) 
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basis function associated with an interior node i. Within the triangles attached to node has 
a vector direction parahel to the edge opposite to node i and forms a loop around node i. 

One the other hand, the tree basis consists of RWG functions that lie along a tree structure 
connecting the centroids of adjacent triangular patches. Fig. S] shows one possible choice of the tree 
basis for this particular triangular mesh. 

















V 








V 











Figure 4: A possible tree basis illustration. Every arrow stands for an RWG function that is the 
member of tree basis. 



3 Proposed method for Neumann problems 

For Neumann problems, the equation in question becomes 

V • e^(r)V(/)(r) = -p(r)/eo for r e ft 

(10) 

£-(/)(r) = ^(r) for r G r 

where F includes all the boundaries of Q. 

A necessary condition for the existence of a solution to this problem is that the source and boundary 
data satisfy the compatibility condition [30] 

/ -e{r)g{r)dl^ [ p{r)dr = 0. (11) 
Jr Jn 

For Neumann problems in homogeneous media, the basic idea of proposed method and the solution 
process have been reported in [T9l. In this section, we extend this approach to handle more complex 
Poisson problems that involve inhomogeneous media. 

3.1 Algorithm outline 

We outline the overall solution procedures as follows: 

1. Acquire the tree space part of electric flux density: By solving a matrix system resulting 
from V • D = p, we first obtain the tree space part of electric flux density, D^, that corresponds to 
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the second term on the right side of Eq. (j8j). Taking advantage of the fast tree solver, we can obtain 
Y>t in 0{Nt) operations. 

2. Field projection onto loop space: In order to get the first term on the right side of Eq. (|8]), 
the loop-space part D/, we need to follow up with a projection procedure. Since the electric field E is 
curl free and orthogonal to the divergence-free space, this can be achieved by projecting the electric 
field E onto the loop space and orthogonalizing it with respect to the loop space, as presented in 
Section [331 By this technique, the desired electric field will be retrieved. 

3. Solution of potential: The potential distribution can be found by solving another matrix 
system accelerated by the fast tree solver that is of 0(7Vt), in a manner similar to that of the first 
procedure. 

3.2 Solution of 

We start with the loop-tree decomposition of electric flux 

Ni Nt 

D(r) = ^/,L,(r) + ^t,T,(r) (12) 

i=l i=l 

where the first term of right hand side is the loop-space part, and the second one is the tree space 
part. 

Meanwhile, from electrostatic theory, we have 

V • D(r) = p{r). (13) 
The charge density /)(r) can be expanded in terms of pulse expansion functions namely, 

p{r) = J2QnPn{r) (14) 



where Np is the number of triangle patches. 

When Eqs. ([12]) and ([Tlj) are substituted into Eq. (p!3|) and testing Eq. ([13]) with a set of pulse 
functions, we have a matrix system 

K-It=Vp (15) 

where 

[K],. = (Pi(r),V-T,-(r)) (16) 

[V,l = {Pi{r),p{r))=q, (17) 

and It = [ti t2 ts • • • ^atJ"^. The number of tree basis functions, Nt^ is equal to Np — 1. There is one 
coefficient in ([T4|) that is not needed due to the charge neutrality. In other words, the last coefficient 
can be derived from the front Np — 1 coefficients, i.e., t^p = — Xli^i ^ '^i- 
In the above, the inner product between two functions is defined as (31] 

(/i,/2)= j /i(r)/2(r)dr (18) 

where the integral is assumed to converge. Moreover, the above matrix system is irrelevant to the 
loop-space part of electric flux because it is divergence free, i.e., V • D/ = 0. 
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3.2.1 Impose Neumann boundary condition 



The above derivation only works for zero Neumann boundary condition, namely, 

d 

--0(r) =0 for r G r. 



When the normal components of the electric field is not trivial on the boundary, we need to handle it 
carefully 

Physically, Eq. (pTj) implies that Neumann boundary condition corresponds to the 2D surface 
charge density introduced into the solution system. Hence, we can add extra half RWG basis functions 
at the boundary where non-zero Neumann boundary condition appears. The coefficients of those half 
RWG functions are proportional to normal electric field at corresponding boundary points. 

3.2.2 Fast tree solver 

The matrix system (p!5]) can be solved with 0{Nt) operations using the fast tree solver. Its key 
component relies on back substitution along a tree structure (e.g., see Fig. [4]). From a mathematical 
point of view, we can invert the matrix K of Eq. ([15]) with linear complexity since a row permutation 
of the matrix is a lower triangle matrix with very few elements per row. This can be achieved with 
the help of topological information. We refer interested readers to Chapter 5 in [28 for details. 

3.3 Loop space projection 

If the media is homogeneous, both the electric fiux and electric field have no loop-space or divergence- 
free component existing, namely, VxD = VxE = 0. Hence, we can project the obtained onto the 
loop space or divergence-free (solenoidal) space in order to remove the divergence free part to retrieve 
the desired field. This is because the desired electric field is orthogonal to the loop (divergence-free) 
space. The procedure is called divergence- free field removal as in [19]. 

Unfortunately, it is not the case in inhomogeneous media where the electric fiux has both divergence- 
free and curl- free parts. It is not orthogonal to the loop space any more. However, the electric field 
E is curl free and orthogonal to the divergence- free space. Hence, we project the electric field E onto 
the divergence- free space, in order to remove its divergence- free component. 

By solving the Eq. ([I5j), coefficients {U}^^ have been obtained and known while {k}^^ are the 
unknowns. We can transform Eq. ([12]) to electric field 



Since the electric field lives in the curl free (irrotational) space, which is orthogonal to the loop space, 
we have 



Testing Eq. (fT9|) by loop basis functions and substituting the above into it, we obtain a matrix system 




(19) 



(Li,E)=0. 



(20) 



(21) 
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where 



V, = -R.I, (23) 
[^,, = /L,(r),5#\ (24) 



= [^1, hi • • • 7 ^NiV is unknown, and 1^ = [/i, ^2, ^3, • • • , ^NtV known. 
Many iterative solvers could be chosen to solve Eq. ([21]), such as CG, Bi-CGSTAB [32j and GM- 
RES [33l[3l]. According to our experience, GMRES has the best performance if no preconditioning 
technique is applied. Once 1/ is found, the numerically approximated curl-free E field is known. From 
it, we can derive the potential (j). 

3.4 Solution of potential 

According to the electrostatic theory, the potential satisfies 

- V(A(r) = E(r). (25) 
We can further expand it by a set of pulse functions 

(/)(r) = ^^,P,(r). (26) 

By substituting Eqs. (p!9|) and ([26|) into Eq. (|25]) and testing it with a set of pulse basis functions, 
it gives 

• = (27) 



where 



[V4 = -(T.(r).^') (28) 



" — T 

K 



= {P,(r),V-T,(r))= [KJ^.^. (29) 

and = [vi 1^2 1^3 ' ' ' i^Np-i]^ ' In the same manner as Section [321 there is one coefficient in (|26|) that 
is dependent on the others. Such coefficient has to be specified as a particular value before Eq. (|27|) 
is solved. This is because the solution of the Neumann problem is found up to an arbitrary constant, 
which corresponds to a reference potential physically. This reference potential does not affect the 
solution of the pertinent problem. This is also reffective of the fact that the inverse of the gradient 
operator is not unique. 

It is interesting that the resulting matrix is just a transpose of the matrix K from Section 13.21 
Hence, we can solve (|27|) with the fast tree solver in 0{Nt) operations similar to Section [321 
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4 Special treatment for Dirichlet boundary 

Handling Dirichlet boundary conditions is of importance since Poisson problems with Dirichlet bound- 
ary conditions commonly occur in practical applications. As discussed in above section, our proposed 
Poisson solver works for Neumann problems. When it comes to the Dirichlet problems or mixed 
boundary problems, however, we need to transform the original problem to the Neumann problem to 
approximate it numerically. 

If we consider the Poisson problem as shown in Fig. [H then the governing equation is 

V • e^(r)V^(r) = -p(r)/eo for r e ft 

I 

(30) 



£(^(r) = ^(r) for r G Fat 



(/>(r) = Vi for r G Tdi 

^(r) = Vr for r eTor 

where Vi and Vr are potential values imposed on left and right part of Dirichlet boundary, F di and 
Tor^ respectively. 

At the first place, we transform the above to be the superposition of two Poisson problems with 
appropriate boundary conditions. The first problem (denoted by PI) is 

V • e^(r)V^i(r) = -/>(r)/eo for r e Q 
^Mr)=9{r) forrGF^ 

(31) 

^i(r) = Vi for re Tdi 

M^) = V; for re Tor 

where is a potential value arisen when we approximate this set equations. 
The second one (denoted by P2) is as follows 

V-e^(r)V^2(r) =0 for r G 
£(/)2(r)=0 forrGF^ 

(32) 

(j)2{r) =0 for r G F^^^ 

M^) = Vr-V; for r GFi,,. 

Obviously, the solution of original problem is just 

(/)(r)=(/)i(r) + (/)2(r). (33) 

Although these two sets of equations satisfy the Dirichlet boundary condition, we can resolve them 
by finding the solutions of two Neumann problems, as described in the following. 
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4.1 Solution of PI 




Figure 5: A schema for treatment of Dirichlet boundary conditions: Two extended regions, and 
Qr, are introduced and assumed filled with extremely high permittivity material. 

From electromagnetic theory, the surface of a perfect electric conductor (PEC) or perfect magnetic 
conductor (PMC) serves as equipotential surface under static or quasistatic conditions. Hence, a 
Dirichlet boundary condition occurs when the equipotential value is known. In addition, it is known 
that a PEC may be regarded as a dielectric material with infinite permittivity. 

Based on the above theory, a technique has been developed to surmount the difficulties arising from 
Dirichlet boundary conditions. Its basic scheme is illustrated in Fig. O Two small regions, namely, 
the left extended region Ql and the right extended one Qr^ are introduced into the solution system; 
both regions are filled with enormously high permittivity dielectric materials that mimic PECs to 
achieve equipotential surfaces. In the context of numerical implementation, we normally endow these 
materials with an enormously large relative permittivity, e.g., 10^. 

Since the electric field inside high permittivity material becomes weak, and the electric field is 
related to the gradient of potential (E = — V0), we can reasonably assume that the boundary of 
extended regions, F^, has homogeneous Neumann boundary condition, where 

F'^ = {d^L - Tdi) + {d^R - Tor) 

as illustrated in Fig. O 

As a result, we obtain the following equation subject only to the Neumann boundary condition to 
resolve problem PI: 

V • e^(r)V^i(r) = -pi(r)/eo for r e Q 

£0i(r)=O forrGFi, (34) 

^(/)i(r) = ^(r) for r G Fat. 

It should be mentioned that, the charge distribution pi(r) has identical distribution as the original 
p(r) in 1], whereas the extended regions, i.e., and Qr, need to be treated carefully as below. 
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At first, in accordance to (fTTj) , pi(r) must satisfy 

/ pi(r)dr= / e{r)g{r)dl- [ p(v)dv = Q' , (35) 

wiiere Q' is just a constant for a particular problem. It implies the quantity of charge for which the 
solution should compensate to guarantee charge neutrality. To this end, a simple but effective way is 
to allocate a line charge with the quantity of Q' inside either 0.^ or 

Next, with appropriate setting of pi(r) in extended regions, we can solve Eq. by our new 
method in Section [3l It should be noted that, when solving the above equations, we must adjust 
the reference potential so that ^i(r) has the value of Vi on the left extended region VLl- This can 
be done because the Neumann problem needs a specified a reference potential point, as discussed 
in Section 13. 4[ and high permittivity material yields an equipotential on the surface. Hence, we can 
adjust that reference potential for convenience. In the right extended region VLr^ by contrast, a floating 
potential F/, that is not necessarily equal to F^, will appear (we call it floating potential because this 
equipotential value is unknown before the solution is obtained). The arising of this floating potential 
does not matter for our method since we can compensate it in the following procedure. 

4.2 Solution of P2 




Figure 6: A schema for solving P2: Two extended regions, VLl and VLr^ are introduced and assumed 
filled with extremely high permittivity material; two line sources that have identical quantity but 
opposite signs are allocated vhVLl and VLr respectively. 

We have shown in the above section that the Dirichlet problem PI can be solved by a Neumann 
problem with extended high permittivity regions. This technique can be applied to solving the second 
problem P2 as well. The solution of P2 is obtained by solving 

V . e,(r)V(/)2(r) = Q[5(v - n) - 5(v - V2)] r G (1^ + (^l + ^r) 

(36) 

|^(/>2(r)=0 forrGFi^orF^ 

where 5{y) refers to a 2-dimensional Dirac delta function that corresponds to a line charge and Q is 
a value relative to the potential difference Vr — V^. between VLl and VLr. Moreover, ri is one point in 
Q.L while Y2 is another point of VLr. This scheme is illustrated in Fig. [6l 
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A numerical difficulty, however, arises in the above equations because the value of Q is not deter- 
mined at this point. Fortunately, this value of Q is unimportant and can be ignored in the solution 
process. To this end, instead of solving Eq. (|36|) directly, we first find the solution of another Poisson 
problem with unit line charge: 

V • e^(r)V^2(r) = 6{r - n) - 6{r - rs), r G (1^ + I^l + ^r) 

(37) 

£(^2(r) =0, for r gT^ or T^. 

Next, having obtained the solution ^2, we can examine the potential difference between two ex- 
tended regions, and Qr. More specifically, suppose this potential difference is V2, we have 

Mr)=^^^^Mr)- (38) 

This is because for Eqs. (|36|) and (|37|) the potential difference between two ends are linearly propor- 
tional to the quantity of charge in question. 

Clearly, the solution process of Eq. (|37j) has no difference from that of PI. In addition, since 
we need to solve this equation only once for a particular geometry, the computational load of this 
procedure is small in many applications, such as the electron transport problem in which Poisson 
equations are solved repeatedly. 

5 Numerical results 

In this section, several examples of 2-D Poisson problems will be solved to demonstrate the effectiveness 
of the proposed method. First, two typical 2-D Poisson problems will be simulated to show the validity 
and efficiency. Then, this solver will be applied to a Poisson problem that arises from a practical 
double-gate MOSFET electron transport simulation. 

This new algorithm has been implemented in C++ and compiled with an Intel compiler. Moreover, 
all simulations listed below are performed on an ordinary PC with 2.66 GHz CPU, 4 GB memory and 
Windows operating system. Furthermore, in order to give a fair comparison, the same sparse matrix 
structure is used for both the proposed Poisson solver (denote by PPS below) and the finite element 
method (FEM). 

5.1 Simple heterogeneous Poisson problem 

First, we consider a 2-D Poisson equation in heterogeneous media as a reference problem: 

V • er(x, ?/)V^(x, ?/) = — 7rcos(7rx) — 7rcos(7r?/) {x,y)eQ (39) 

where ft = [0, 1] x [0, 1], and the relative permittivity 




1, x<0.5 

2, x>0.5. 
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For convenience, we assume that this problem has zero Neumann boundary condition. With a 

reference potential 2/7r imposed at the origin, it has the close form solution 

. , . cosfyrx) + cosfyr?/) 
<\>\^.V) = 7 ^ 




Figure 7: The electric electric flux density calculated by the proposed method. Left: x component of 
the electric flux density D. Right: y component of the electric flux density D. 
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Figure 8: Comparison of the electric flux, D, along one line between result of proposed Poisson solver 
(PPS) and that of FEM. Left: x component of the electric flux. Right: y component of the electric 
flux. 

Fig.[7]shows calculated x- and ^/-components of electric flux density D. Because there is discontinu- 
ity for at X = 0.5, ^/-components, Dy^ appears as an abrupt change in the middle correspondingly, 
which is in complete agreement with the fundamental theory of electromagnetism. Therefore, our 
proposed method has the capability to handle discontinuous media. In addition, further comparisons 
between electric fluxes obtained from proposed method and traditional flnite element method are 
given in Fig. [H from which we can see a good agreement at every sample points. Moreover, Fig. [9] 
shows calculated potential distributions, in which the result of proposed method is given in the left 
flgure while the FEM one is given in the right flgure as a reference. 

To examine the computational complexity, we plot the total computing time as the mesh density 
increases (see Fig. [10]). This computing time includes three parts: (a) the flrst fast-tree solution 
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Potential distribution 



Potential clistribution(FEM) 





Figure 9: Calculated potential distribution. Left: the result of the proposed method. Right: the 
result of FEM. 



Computational complexity at 0.01 stoping criterion 



X PPS-GIV1RES(60) 

0{N) 

V FEM-GMRES(60) 
0(N^-^) 



Number of triangles 



Computational complexity at 0.001 stoping criterion 



X PPS-GIVIRES{60) 

0(N) 

V FEM-GMRES{60) 
0(N^-^) 



Number of triangles 



Figure 10: Computational complexity comparison between the proposed Poisson solver (PPS) and 
FEM method. Left: 0.01 stopping criterion. Right: 0.001 stopping criterion. 



time corresponding to Section 13. 2| (b) loop-space projection time corresponding to Section 13.31 and 
(c) the second fast tree solution time corresponding to Section 13. 4[ Among them, (a) and (c) can 
be computed expeditiously with veracious 0{N) complexity which has been proved in [25 , whereas 
(b) includes an iterative procedure that dominates the total solution time. Apparently, the iteration 
solution time depends on the iterative solver type and the required accuracy level. In our numerical 
experiment, we observe that the total solution time of our proposed method is close to 0{N) complexity 
when GMRES solvers, without any preconditioning techniques, are applied in this problem with two 
different stopping criterion: 1 x 10~^ and 1 x 10~^. By contrast, as shown in plots of Fig. [TOl the 
computational complexities of traditional FEM method are approximately of 0(7V^-^) and 0{N^'^) at 
the given stopping criterion, respectively. Hence, this new proposed method demonstrates impressive 
efficiency compared with the FEM method. 
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(-0.7,1.7) 




(-2,0) (2,0) 



Figure 11: A two dimensional region where the Poisson problem is defined. 

5.2 Complex Poisson problem with Dirichlet boundary condition 

Next, to show the capability of handling Dirichlet boundary conditions, we simulate a general two- 
dimensional Poisson equation, which is 

V-e^(r)V(^(r) = -J(r-rO for r e ft (40) 

with boundary condition 

^(r) = 1.0 re Fi 
0(r) =0.8 re F2 
-^(f){r) = r G other boundaries, 

and is the point (—0.2,0.6). 

Fig. [TT] shows the specifications of solution region Q. This problem is excited by a line source that 
is located at point with unit charge and imposed Dirichlet boundary conditions on the left and right 
edges. 



Potential distribution Potential distribution(FEM) 




Figure 12: Potential distribution calculated from two method. Left: the proposed method. Right: 
FEM method. 

With the technique described in Section IH we can solve this problem by the proposed method. In 
our numerical experiment, the region is discretized into 133,388 triangle patches, using the following 
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Table 1: A comparison of solution times between proposed method and traditional finite element 
method. 





Proposed method 


FEM 


Stopping criteria 


Iterative steps 


Solution times (second) 


Iterative steps 


Solution times (second) 


1 X 10-2 


26 


0.779 


31 


0.889 


1 X 10-^ 


104 


4.742 


256 


11.216 


1 X 10-"^ 


295 


14.446 


2468 


111.976 


1 X 10-^ 


593 


28.89 


4813 


297.539 



function to approximate the line source 

! 625 \x\ <Om,\y\ <0.02 
otherwise. 

As for the iterative solver, GMRES with restart number 60 is employed. Finally, the calculated 
potential distribution is shown in Fig. [121 The left figure shows the potential distribution calculated 
by our proposed method while the right one gives FEM results as the reference. It is shown that 
our proposed method can produce the same results as the reference one. Hence, not only Neumann 
boundary conditions, but also Dirichlet boundary conditions can be handled by our proposed solver 
competently. 

To inspect the efficiency, we observe the pertinent computing resource of both the proposed method 
and the traditional FEM method. Both methods apply GMRES iterative solver with a restart number 
60, without any preconditioning. Table 15.21 lists the solution time and iteration number for both 
methods. As seen from this table, both iteration number and solution time of the proposed method 
are considerably less than that of the traditional FEM method. 

5.3 Application to double-gate MOSFET 

Finally, the solver is used to simulate a multigate silicon MOSFET, as shown in Fig. [131 which 
is a promising candidate for the next generation nanotransistor. The structure is infinite in the z 
direction. Gate length is denoted by Lg; source and drain extension lengths are denoted by Ls and 
Ld^ respectively; silicon channel thickness is Wch] and oxide thickness is Wqx- In our simulation, the 
device parameters are Lg = 10 nm, Lg = = 4: nm, Wch = 5 nm, and Wqx = 1 nm; doping density 
is = 10^^/m^; the relative permittivity of silicon is 11.9; and relative permittivity of the silicon 
dioxide is 3.8. Other parameters can be found in [7]. 

In this simulation, we find the solution of the following Poisson equation: 

V • [e(x, y)\/VD{x, y)] = q[n{x, y) - Nd{x, y)] (41) 



18 



Ls Lg Ld 




Figure 13: Two-dimensional view of the n-type double-gate silicon MOSFET. 

where Vd is the potential to be found, n(x, y) is electron density distribution given by Schrodinger 
solver. Here, is the doping density and q is the electron charge which is equal to 1.62 x 10~^^ 
coulombs. The simulation domain is the region enclosed by dashed line in Fig. [131 Moreover, the 
Dirichlet boundary condition is enforced at the gate region, whereas the floating boundary condition, 
i.e., the zero normal derivative, is applied at the remaining boundary. 



E 



2 




-8 -6 -4 -2 2 4 6 8 

X(nm) 



Figure 14: Two-dimensional plot of the electron density distribution. 

Fig. [m demonstrates the electron density distribution provided from Schrodinger solver. We then 
use this electron density distribution as the source of Eq. (|4T]) . By our proposed method, the potential 
distribution is the same as the result in [7] and is shown in Fig. [151 

6 Conclusions 

We have proposed and developed an efficient numerical method for solving Poisson problems. It can 
deal with Poisson equation with both Neumann and Dirichlet boundary conditions in a 2-D irregular 
region filled with homogeneous or inhomogeneous dielectrics. In this method, the electric fiux is 
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X(nfin). 



Figure 15: Two-dimensional plot of calculated potential distribution. 

first found as opposed to traditional methods. This can be done by quasi-Helmholtz decomposition 
with loop-tree basis functions whose coefficients could be gained rapidly with fast solvers. Although 
the Dirichlet boundary condition appears as an numerical obstacle for the proposed new method, it 
can be overcome using special treatment: small extended regions with extremely high permittivity 
material are introduced to force the potential to be constant value to imitate the pertinent Dirichlet 
boundary conditions. Through numerical examples, it is shown that this new method can solve a 
general Poisson equation arising from electrostatics robustly, and has better performance than the 
traditional finite element method. According to our numerical experiments, it is observed that the 
computational complexity of this new method is close to 0{N). This method is a novel fast Poisson 
solver that can serves as a feasible alternative to multigrid scheme for Poisson solutions. 
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